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ABSTRACT 


The  two-dimensional  (2D)  spectral  properties,  including  the  degree  of  isotropy,  of 
under-ice  roughness  in  the  Beaufort  Gyre  were  investigated.  Under-ice  roughness  plays 
an  important  role  in  turbulent  ocean  heat,  salt,  and  momentum  fluxes  that  determine  the 
delicate  balance  between  surface  forcing  and  the  ocean  interior  that  sustains  or  reduces 
the  perennial  ice  cover.  Accurate  characterization  of  this  roughness  is  important  for 
numerical  modeling  and  prediction  of  the  Arctic  air-ice-ocean  system,  which  will  play  a 
significant  role  as  the  US  Navy  increases  its  strategic  presence  in  the  ice-diminished 
Arctic  Ocean.  Ice  draft  data  from  moored  upward  looking  echosounders  of  the  Beaufort 
Gyre  Exploration  Project  and  ice  motion  derived  from  trajectories  of  International  Arctic 
Buoy  Program  buoys  were  combined  to  create  directional  sections  of  ice  draft. 
Autocorrelation  functions  were  created  for  each  segment  which  were  then  combined  to 
make  a  2D  autocorrelation  function.  The  2D  FFT  was  taken  of  this  in  order  to  yield  a  2D 
spectrum,  from  which  isotropy  and  other  spectral  properties  could  be  detennined.  The 
processing  routine  was  tested  on  several  synthetic  test  datasets  and  in  each  case  produced 
the  correct  output.  Two  years  of  data  from  three  moorings  in  the  Beaufort  all  provided 
2D  spectra  which  showed  weak  anisotropy  of  the  under-ice  roughness.  Error  estimates 
were  detennined  by  Monte  Carlo  simulations.  For  each  dataset,  most  of  the  variance  of 
the  underside  morphology  was  contained  in  a  wavenumber  band  between  0.001  and 
0.01m'1  (lOO-lOOOm).  There  were  no  distinct,  preferred  wavelengths.  Spectral  levels 
ranged  from  1.2  -  4.2  x  104  m2  for  Mooring  C  (2004-2005)  to  0.2  -  1.6  x  104  nr  for 
Mooring  A  (2004-2005). 
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I.  INTRODUCTION 


A.  CURRENT  EVENTS 

Recent  changes  in  the  Arctic  Ocean  and  specifically  its  diminishing  ice  cover 
have  received  significant  attention  both  in  scientific  circles  and  amongst  the  general 
public.  The  current  International  Polar  Year  is  bringing  increased  awareness  and 
gathering  crucial  data  needed  to  understand  these  changes.  An  ice-diminished  Arctic  is 
also  of  strategic  significance  to  the  national  security  of  the  United  States.  The  viability  of 
both  the  Northern  Sea  Route  and  the  Northwest  Passage  as  navigable  shipping  lanes  for 
the  first  time  in  recorded  history  combined  with  the  new  feasibility  of  resource 
exploitation  (minerals,  oil,  fisheries)  in  the  Arctic  Basin  itself  has  led  many  nations  to 
reassert  or  even  stake  new  claims  in  the  region.  Russia,  Norway,  Sweden,  Finland, 
Denmark,  Iceland,  Canada  and  the  United  States  all  have  major  interests  in  the  Arctic  and 
each  seeks  to  extend  their  Exclusive  Economic  Zone  (EEZ)  in  any  way  possible.  A  130 
year  long  dispute  continues  between  the  United  States  and  Russia  over  the  location  of  a 
line  dividing  possession  of  the  Bering  Sea  between  Siberia  and  Alaska  (Damron,  1997). 
Norway  and  Russia  have  competing  claims  to  fishing  rights  in  the  Barents  Sea  (BBC 
News,  2005).  Canada  insists  that  the  Northwest  Passage  lies  within  Canadian  territorial 
waters  and  that  they  will  assert  sovereign  control  of  it  whereas  the  United  States  and 
many  other  countries  claim  international  rights  to  use  it  for  shipping  (Shukman,  2007). 
Most  recently,  in  the  summer  of  2007,  the  Russian  Federation  submitted  a  formal  claim 
for  over  45%  of  the  Arctic  seafloor,  including  the  North  Pole  itself. 

The  summer  of  2007  also  saw  several  military  developments.  A  Russian 
submersible  symbolically  planted  a  flag  on  the  ocean  floor  at  the  geographic  North  Pole. 
Canada  announced  that  it  would  enhance  its  military  presence  in  the  Arctic,  including 
building  two  new  military  facilities;  one  for  cold-weather  warfare  training  at  Resolute 
Bay  and  the  other  a  deep-water  port  on  the  northern  end  of  Baffin  Island  (Reid,  2007). 
The  Canadian  government  also  placed  an  order  for  six  to  eight  new  Polar  Class  5  Arctic 
Offshore  Patrol  Ships.  The  vessels  will  be  anned  and  according  to  Prime  Minister 
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Harper,  they  will  “allow  us  patrol  the  Arctic  waters  in  the  Northwest  Passage  when  it 
matters  (CBC  News,  2007).”  An  American  guided  missile  cruiser,  USS  NORMANDY 
(CG  60),  conducted  a  patrol  across  the  Arctic  Circle  north  of  Iceland  to  show  that  the  US 
Navy  maintains  the  capability  to  operate  in  the  far  north.  USCGC  HEALY  conducted  a 
bathymetric  survey  of  the  Chukchi  Cap  in  coordination  with  NOAA’s  Office  of  Coast 
Survey  “to  better  understand  its  morphology  and  the  potential  for  including  this  area 
within  the  United  States’  extended  continental  shelf  under  the  United  Nations  Convention 
on  Law  of  the  Sea  (Elliott,  2007).”  In  this  “gold  rush”  atmosphere,  there  is  developing 
potential  for  anned  conflict  over  both  resources  and  access.  At  the  very  least,  the  US 
Navy  can  be  expected  to  conduct  routine  Freedom  of  Navigation  demonstrations,  sending 
ships  to  physically  challenge  territorial  claims  that  the  US  does  not  recognize. 

In  order  to  safely  operate  in  the  Arctic,  whether  a  vessel  is  military,  merchant,  or  a 
tourist  cruise  ship,  ice  modeling  and  the  resulting  forecasting  will  become  as  important  as 
routine  weather  forecasting  is  today.  There  are  a  number  of  challenges  in  ice  modeling. 
Due  to  the  remote  and  harsh  Arctic  environment,  these  difficulties  include  sparseness  of 
in  situ  data.  Furthermore,  there  is  a  limited  amount  of  remotely  sensed  data.  Of  all  of  the 
unknowns,  one  that  is  particularly  challenging  is  morphology  of  the  underside  of  the 
icepack  which  actually  is  in  contact  with  the  ocean.  This  cannot  be  determined  remotely 
with  current  technology. 

B.  UNDER  ICE  ROUGHNESS  FEATURES  AND  DEFINITIONS 

Pressure  ridges  are  fonned  in  the  Arctic  Ocean  when  two  ice  floes  are  forced 
together  by  wind  or  current.  During  the  collision,  areas  of  thin  ice  between  the  floes  are 
crushed  and  sometimes  even  the  edges  of  the  floes  themselves  fracture.  This  material 
accumulates  and  piles  up  so  as  to  form  a  ridge  visible  from  above  the  ice.  An 
approximate  mirror  image  (but  much  larger  in  size)  also  forms  as  an  inverted  ridge 
protruding  downwards  into  the  ocean  away  from  the  underside  of  the  ice  floe.  This 
submerged  feature  is  referred  to  as  an  ice  keel.  It  has  been  estimated  that  pressure  ridges 
and  their  associated  ice  keels  account  for  approximately  half  of  the  ice  volume  in  the 
Arctic  Basin  (Davis  and  Wadhams,  1995). 
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c. 


SIGNIFICANCE 


Ridging  events  influence  the  magnitude  and  direction  of  ice  motion.  Pressure 
ridges  cause  an  increase  in  the  aerodynamic  drag  at  the  air-ice  interface  and  their 
associated  keels  cause  an  increase  in  hydrodynamic  drag  at  the  ice-water  interface.  If  an 
ice  keel  is  deep  enough  it  can  generate  internal  waves  at  the  pycnocline,  exerting  much 
higher  drag  forces.  This  can  lead  to  the  phenomenon  referred  to  by  Ekman  as  “dead 
water”  (Wadhams,  2002).  This  occurs  in  highly  stratified  areas  where  much  of  the 
energy  used  by  a  ship  for  propulsion  is  lost  due  to  the  formation  of  internal  waves  as  the 
keel  plows  though  the  pycnocline.  This  loss  of  momentum  to  internal  wave  generation 
complicates  modeling  of  the  total  air-ice-ocean  momentum  budget,  especially  when  the 
size  and  freqency  of  the  keels  is  unknown. 

It  is  also  important  to  note  the  flow  effects  induced  by  ice  keels.  The  ratio  of  a 
keel  depth  to  the  depth  of  the  oceanic  boundary  layer  is  very  large  and  can  even  exceed 
one,  as  was  the  case  at  the  Applied  Physics  Laboratory  Ice  Station  (APLIS)-2007.  This 
makes  a  field  of  large  ice  keels  analagous  to  a  large  mountain  range  in  the  atmospheric 
boundary  layer.  Keels  can  cause  flow  acceleration,  blocking  and  stagnation  along  with 
turbulent  fluxes  of  heat,  salt,  and  momenturm.  Simulations  indicate  that  the  turbulent 
wakes  generated  as  ice  keels  are  pushed  through  the  water  can  exend  several  hundred 
meters  downstream.  Large  Eddy  Simulation  modeling  has  shown  that  a  10m  deep  keel 
can  cause  a  five-fold  increase  in  heat  flux  (Skillingstad  et  al,  2003). 

There  are  engineering  applications  which  need  to  account  for  ridging.  Ice  ridges 
are  the  deciding  factor  in  determining  strength  requirements  for  Arctic  offshore  platforms 
and  icebreaker  hull  design.  The  depth  of  ice  keels  and  their  potential  to  scour  the  bottom 
in  shallow  water  determines  the  depth  of  pipeline  burial  (Wadhams,  2002). 

There  are  significant  naval  applications  as  well.  Safety  of  submarine  navigation 
in  the  shallow  waters  of  the  Arctic  shelf  and  ability  of  a  submarine  to  surface  through  the 
ice  are  affected  by  ice  keels.  Possibly  even  more  important  are  the  acoustic  effects 
attributed  to  ice  keels.  Due  to  its  density  structure,  the  Arctic  Ocean  is  characterized  as 
an  upward  refracting  acoustic  environment  (Fig.  1).  Acoustic  energy  gets  turned  towards 


3 


the  ice  where  it  hits  and  is  scattered  as  a  function  of  the  ice  bottom  roughness  (Fig.  2). 
An  even  more  extreme  case  is  where  an  ice  keel  is  deep  enough  to  cause  blocking  of  the 
acoustic  path  (Bourke  and  McLaren,  1992).  The  acoustic  shadows  cast  by  an  extremely 
deep  keel  make  it  possible  for  a  submarine  to  hide  in  its  lee,  shielding  it  from  active 
sonar.  Also  of  significance  to  submarine  operations  is  the  ambient  noise  of  the  ice  pack 
itself.  High  levels  of  ambient  noise  are  associated  with  ice  floe  collisions,  especially  at 
the  ice  edge  (Wadhams,  2002).  This  phenomenon  could  greatly  complicate  passive 
sonar  detection. 

SOUND  VELOCITY  (M/ SEC) 


1430  1470  ISO  1600 


Figure  1 .  Simplified  Arctic  sound  speed  profiles.  Model  1  is  based  on  a  sequence  of 
plane  parallel  layers  each  with  a  constant  gradient  whereas  Model  2  is  based  on  a  series 
of  flat  layers  of  constant  velocity  and  density  (Kutschale,  1969). 
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Figure  2.  Acoustic  ray  paths  for  a  source  at  100m  depth  using  Model  1  from  Fig.  1 
(Kutschale,  1969). 

The  effects  of  ridging  may  extend  to  mechanisms  important  to  climate  change. 
The  more  ridged  the  ice  is,  the  greater  the  volume  of  fresh  water  that  can  cross  flux  gates 
such  as  Fram  Strait.  Ridging  also  increases  the  melting  rate  of  ice,  by  increasing  the 
surface  area  of  the  ice  in  contact  with  the  wanner  ocean  underneath.  If  an  ice  keel  is 
deep  enough,  it  can  cause  deeper  (wanner)  water  to  be  dragged  up  behind  it  and 
subsequent  turbulent  mixing  will  bring  this  warm  water  into  contact  with  the  underside  of 
the  ice,  further  melting  the  ice  from  below.  Observational  data  support  the  assertion  that 
ridged  ice  melts  faster  than  undefonned  ice  (Davis  and  Wadhams,  1995). 

D.  STUDY  OBJECTIVES 

The  purpose  of  this  study  was  to  produce  a  two-dimensional  (2D)  spectral 
respresentation  of  ridge  properties  and  roughess  of  the  underside  of  Arctic  sea  ice  which 
would  allow  determination  of  anisotropy.  The  spectral  approach  makes  it  possible  to 
determine  preferred  direction,  wavelength,  and  magnitude  of  the  roughness.  If  the 
roughness  is  isotropic,  coupled  air-ice-ocean  models  can  reliably  use  a  single  coefficient 
for  roughness,  regardless  of  the  direction  of  the  wind-forced  ice  motion.  However,  if  the 
roughness  is  anisotropic,  directionally  variable,  direction-dependent  coefficients  may  be 
necessary  to  accurately  simulate  ocean-ice  exchange  rates.  This  study  combines  moored 
upward  looking  sonar  (ULS)  derived  ice  draft  data  with  motion  determined  from  drifting 
buoys  to  sparsely  sample  and  develop  2D  spectra  of  under-ice  roughness  topography. 
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II.  DATA  SELECTION 


A.  ICE  DRAFT  DETERMINATION 

Significant  efforts  have  been  made  to  infer  under-ice  topography  from  airborne 
laser  techniques.  There  have  been  simultaneous  submarine  upward  looking  sonar  and 
and  LIDAR  equiped  aircraft  transits  that  attempted  to  establish  relationships  between  the 
two  surfaces.  While  there  was  some  success,  the  amount  of  snow  cover  is  variable  and 
difficult  to  determine,  thus  contaminating  the  freeboard  to  draft  ratio  for  the  ice  alone 
(Wadhams,  2002).  A  typical  ratio  used  as  a  “rule  of  thumb”  for  ice  keel  draft  to  ridge 
height  is  3.2:1.  Recently,  measurements  of  a  ridge  at  APLIS  2007  did  not  follow  the 
relationship.  A  hole  was  melted  through  the  ice  near  an  ice  ridge  to  deploy  instruments. 
The  ice  was  approximately  13.5m  thick,  which  suggests  that  the  hole  went  through  the 
ice  keel  itself.  The  associated  pressure  ridge  had  a  maximum  surveyed  height  of  2.88m 
above  sea  level,  implying  that  an  its  underlying  keel  would  be  on  the  order  of  9.2m  depth 
(Fig.  3).  Even  if  the  tip  of  the  keel  had  been  offset  from  the  ridge  through  some  sort  of 
shearing  or  rafting  event,  it  still  would  not  closely  follow  the  ratio.  Furthermore,  any 
mirror  image  of  ridge  and  keel  geometry  would  to  some  extent  be  distorted  for  multiyear 
ice,  which  has  undergone  annual  cycles  of  melting  and  freezing  cycles,  along  with 
mechanical  stress  and  strain  events  and  their  subsequent  deformation.  This  observation 
shows  significant  differences  between  ice  keel  and  surface  ridge  structure. 
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Figure  3.  Surveyed  ice  ridge  from  APLIS  2007  and  an  extrapolated  ice  keel  using 
the  stadard  ratio  for  ice  keel  draft  to  ridge  height  of  3.2: 1 . 

There  are  two  common  methods  of  obtaining  reliable  ice  draft  distribution  data, 
both  of  which  use  ULS.  One  method  is  to  mount  an  ULS  on  a  submarine  which  then 
operates  under  the  ice,  recording  ice  draft  as  it  goes.  The  second  method  is  to  mount  the 
ULS  to  a  bottom  mooring  which  extends  to  within  about  50m  of  the  surface  and  record 
the  ice  draft  which  comes  into  its  field  of  view.  There  are  strengths  and  weaknesses  in 
both  methods. 

1.  Submarine  Data 

The  strength  of  submarine  data  is  that  it  can  give  a  near-instantaneous  snapshot  of 
a  broad  area  because  the  submarine  can  move  quickly  in  relation  to  the  speed  of  ice 
motion.  During  submarine  missions,  ULS  data  storage  and  power  requirements  are  not 
major  concerns. 
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However,  submarine  data  have  a  number  of  limitations.  Data  can  be 
contaminated  by  vessel  depth  or  speed  variations.  Submarine  cruises  are  spatially  and 
temporally  sparse.  The  areas  surveyed  are  rarely  if  ever  repeated.  Furthermore,  except 
for  the  SCICEX  missions,  collecting  environmental  data  is  not  the  submarine’s  primary 
mission  and  therefore  quality  may  suffer  somewhat. 

Submarine  ULS  data  are  rarely  released  publicly  due  to  security  issues  concerning 
sovereignty  of  territorial  waters.  Submarine  datasets  were  rarely  released  to  the  public 
before  A1  Gore,  then  Vice  President,  petitioned  to  declassify  data  in  an  area  which  has 
become  known  as  the  “Gore  Box”  (National  Science  Foundation,  1998).  This  release 
area  has  since  been  expanded  by  the  Chief  of  Naval  Operations  beyond  the  original 
“Gore  Box”  (Fig.  4)  (National  Snow  and  Ice  Data  Center,  2006). 


Figure  4.  Original  “Gore  Box”  and  current  release  data  for  some  submarine 
collected  upward  looking  sonar  ice  draft  data  (National  Snow  and  Ice  Data  Center,  2006). 


There  are  additional  safety  issues  to  consider.  It  is  dangerous  for  a  submarine  to 
be  in  shallow  water  under  the  ice  with  little  room  to  maneuver.  For  example,  USS 


9 


SARGO  (SSN  583)  damaged  her  sail  in  a  collision  with  an  ice  keel  in  the  shallow  and 
confined  waters  of  the  Bering  Strait  during  the  winter  of  1960  (Fig.  5)  (Nicholson,  2008). 


Figure  5.  Crewmembers  inspect  sail  damage  on  USS  SARGO  (SSN-583)  after 
colliding  with  an  ice  keel  in  the  Bering  Sea  in  February  1960  (Nicholson,  2008). 

2.  Mooring  Data 

Mooring  mounted  ULS  allows  for  the  collection  of  a  long  term  record  of  sea  ice 
draft.  The  locations  are  usually  re-occupied  so  multiple  years  of  data  can  be  collected  in 
this  manner.  It  is  cheaper  and  safer  to  deploy  a  mooring  than  to  operate  a  nuclear 
submarine.  Moorings  can  be  located  in  the  very  shallow  waters  of  the  continental  shelf 
where  submarines  cannot  operate.  Finally,  these  types  of  datasets  are  not  classified. 

Moored  ULS  has  weaknesses  too.  Temporal  coverage  is  good  but  spatial 
coverage  is  limited  to  a  single  point  for  each  mooring.  The  unit  must  be  completely  self- 
contained  and  rugged  enough  to  survive  the  harsh  conditions,  placing  restrictions  on 
memory  capacity  and  power  requirements.  Another  weakness  is  that  the  data  cannot  be 
examined  during  the  period  of  collection.  The  failure  of  a  unit  will  not  be  apparent  until 
recovery,  which  is  typically  a  year  after  deployment. 
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After  comparing  the  available  submarine-mounted  and  moored  ULS  data  sets,  it 
was  decided  to  use  the  moored  ULS  data  because  the  observations  were  recent,  year- 
round,  and  available  for  a  number  of  different  locations.  These  attributes  make  seasonal, 
annual,  and  regional  comparisons  of  under-ice  roughness  possible.  Specifically,  the  data 
used  in  this  project  were  from  the  Woods  Hole  Oceanographic  Institute’s  Beaufort  Gyre 
Exploration  Project  (BGEP)  (Proshutinsky,  2008). 

B.  CONVERTING  TO  A  SPATIAL  SERIES  OF  ICE  DRAFT 

The  BGEP  moored  ULS  provide  time  series  of  ice  draft,  z(t).  However, 
determination  of  isotropy  requires  a  2D  spatial  determination  of  ice  draft,  z(x,y).  To 
convert  from  one  to  the  other,  the  time  series  is  multiplied  by  the  sampling  rate  and 
advection  speed  of  the  ice.  Therefore,  the  motion  of  the  ice  needed  to  be  determined. 

C.  ICE  MOTION  DETERMINATION 

Ice  motion  can  be  derived  from  a  variety  of  methods.  These  include  Eulerian 
measurements  with  an  acoustic  Doppler  current  profiler  (ADCP)  pointed  upwards  from  a 
fixed  location,  Lagrangian  motion  detennined  from  drifting  buoys,  and  satellite  remote 
sensing  products. 

Ideally,  there  would  be  an  ADCP  collocated  with  every  ULS  to  give  high 
resolution  local  ice  motion.  Unfortunately,  this  type  of  data  is  not  always  available.  For 
this  project,  ice  motion  was  determined  from  the  position  fixes  of  buoys  coordinated  by 
the  International  Arctic  Buoy  Program  (IABP)  (Rigor,  2002).  This  buoy-derived  motion 
was  preferred  over  a  number  of  satellite  products  such  as  those  derived  from  AVHRR, 
scatterometry,  and  other  merged  datasets  because  of  its  higher  temporal  and  spatial 
resolution  (Fowler,  2003).  Most  satellite  products  give  averaged  12  or  24  hour  motion 
fields.  However,  in  order  to  measure  2D  under  ice  roughness,  the  motion  vectors  had  to 
be  available  with  a  temporal  resolution  on  the  order  of  one  hour  so  that  the  frequently 
dominant  ice  motion  at  these  latitudes,  inertial  oscillations,  could  be  resolved.  This  high 
resolution  data  could  only  be  obtained  from  the  IABP  buoys. 
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III.  DATA  PROCESSING 


A.  ICE  DRAFT  DATA  PROCESSING 

As  part  of  BGEP,  three  ULS  units  were  mounted  on  fixed  moorings  (referred  to 
as  ‘A’,  ‘B’,  and  ‘C’)  under  the  ice  of  the  Beaufort  Gyre  (Fig.  6). 


Mooring  Locations  for  BGEP 
2003  -  2005 
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Figure  6.  BGEP  Mooring  Focations. 

The  moorings  were  anchored  in  water  over  3500  m  deep,  with  the  uppermost 
instrument  on  the  mooring  being  the  UFS,  which  was  located  approximately  50m  below 
the  ice  bottom  (Fig.  7).  The  UFS  instruments  used  were  ASF  Environmental  Sciences 
model  IPSA,  which  operated  at  a  frequency  of  420  kHz  with  a  beamwidth  of  1.8  degrees, 
yielding  a  footprint  on  the  underside  of  the  ice  of  about  2m.  Range  between  the  UFS  and 
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the  ice  bottom  was  sampled  at  2  seconds.  Total  hydrostatic  pressure  and  water 
temperature  were  sampled  every  40  seconds  (Krishfield  and  Proshutinsky,  2006). 


Palm 


Figure  7.  Schematic  of  ULS  mooring  (Krishfield  and  Proshutinsky,  2006). 

In  general,  ULS  data  require  significant  processing  and  despiking  before  they  are 
usable.  For  this  project  in  particular,  ice  draft  was  calculated  by  correcting  the  measured 
range  for  tilt  and  sound  speed  variations  and  then  correcting  the  total  hydrostatic  pressure 
measured  at  the  transducer  for  variations  in  atmospheric  pressure  and  changes  in  seawater 
density.  Sound  speed  and  density  are  only  determined  from  CTD  casts  at  the  beginning 
and  end  of  the  deployment.  During  this  deployment,  it  was  impractical  to  have  a  CTD 
continuously  profiling  so  other  methods  are  required.  Temperature  measured  at  the  ULS 
can  locally  correct  sound  speed  and  density,  while  corrections  can  also  be  determined  by 
using  the  range  observations  during  periods  of  open  water  (where  draft=0 ).  However, 
this  method  has  limitations  because  it  is  often  difficult  to  distinguish  frazil  and  thin  ice 
from  open  water.  Since  in  situ  observations  of  atmospheric  pressure  are  very  sparse  in 
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the  polar  regions,  they  were  obtained  from  NCEP  sea  level  pressure  reanalysis  data. 
After  all  corrections  have  been  made,  draft  calculations  have  an  error  of  ±5 -10  cm. 
(Krish Held  and  Proshutinsky,  2006). 

A  large  source  of  error  in  all  ULS  observations  is  the  finite  beam  width  of  the 
transducer.  Although  it  is  assumed  that  the  draft  value  assigned  to  a  particular  location 
on  the  ice  is  along  the  transducer’s  axis,  this  is  not  always  the  case.  On  sloping  features, 
the  first  return  of  the  sonar  pulse  may  be  due  to  an  off  axis  feature  which  is  shallower. 
This  will  cause  the  apparent  draft  at  a  point  to  be  greater  (Fig.  8).  The  net  effect  is  to 
overestimate  ice  draft  and  to  broaden  the  apparent  width  of  features  (Melling  et  al,  1995). 
The  only  way  to  mitigate  this  effect  is  to  use  an  ULS  with  as  narrow  a  beam  width  as 
possible.  The  IPS-4  transducer  has  a  beamwidth  of  1.8  degrees  which  is  narrow  in 
comparison  to  sonars  used  to  generate  previous  datasets,  which  used  beamwidths  on  the 
order  of  14  degrees  (Wadhams,  1978). 


Figure  8.  Effect  of  echosounder  beam  width  on  apparent  keel  profile.  Note  that  the 
apparent  width  and  average  draft  increase  (Wadhams,  1978). 

Other  sources  of  error  are  due  to  phenomena  occurring  in  leads.  Waves  can  give 
an  apparent  draft  if  sampling  occurs  while  its  trough  is  over  the  ULS.  Likewise,  the  crest 
of  a  wave  can  give  a  solution  with  a  negative  draft,  which  would  physically  mean  the 
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sonar  was  looking  out  of  the  water.  Another  effect  is  from  wave  breaking,  which  can 
inject  clouds  of  bubbles  into  the  water.  These  in  turn  may  appear  as  legitimate  ice  drafts 
in  the  sonar  record  (WCRP,  2004). 

Another  major  complication  in  interpreting  ULS  observations  is  the  effect  of 
pressure  ridge  orientation  relative  to  the  direction  of  sampling.  If  a  keel  is  sampled 
perpendicular  to  its  axis  of  lineation,  there  is  no  problem.  However,  if  the  keel  is 
sampled  at  some  other  angle,  the  effect  is  that  the  keel  appears  to  be  broadened  and  the 
spacing  between  features  of  similar  alignment  is  lengthened  as  well  (Fig.  9  and  10).  The 
wavelength  will  always  appear  lengthened;  it  cannot  appear  shortened.  This  leads  to  a 
“reddening”  of  the  spectra  whereby  actual  wavenumbers  are  shifted  to  lower 
wavenumbers.  The  most  extreme  case  occurs  when  sampling  along  the  axis  of  a  keel, 
where  the  wavelength  or  spacing  between  keels  becomes  infinite  and  wave  number  goes 
to  zero. 
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Ice  Keel  with  Two  Sampling  Directions 
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Figure  9.  (a)  View  from  below  ice  keel  showing  sampling  directions  perpendicular 

to  the  keel  and  cutting  across  it  at  75°.  (b)  Apparent  keel  profdes  as  a  function  of 
sampling  angle.  6  =  0°  and  6  =  90°  correspond  to  sampling  directions  perpendicular  and 
parallel  to  the  keel  axis,  respectively.  Note  that  the  apparent  width  and  average  draft 
increase  as  the  angle  becomes  greater. 
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Field  of  Ice  Keels  with  Two  Sampling  Directions 
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Figure  10.  (a)  View  from  below  field  of  keels  with  sampling  directions  0  and  75°.  (b) 

Apparent  draft  profiles  from  field  of  regularly  spaced  keels  as  a  function  of  sampling 
angle.  Note  the  increase  in  apparent  spacing  of  the  ice  keels  as  the  angle  increases. 
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This  problem,  sometimes  refered  to  as  spatial  aliasing,  needs  to  be  resolved  in 
order  to  generate  a  meaningful  2D  spectrum  of  under-ice  roughness.  As  the  examples 
(Fig.  9  and  10)  above  have  illustrated,  a  ID  slice  through  a  2D  field  is  reddened  (except 
for  the  one  case  where  the  slice  happens  to  be  exactly  perpendicular  to  the  orientation  of 
the  field  of  features).  This  reddening  occurs  because  the  2D  wavenumber  information  is 
lost.  If  the  whole  2D  spatial  field  were  available,  obtaining  a  correct  2D  spectrum  would 
be  easy  with  a  2D  FFT.  This  project  used  a  processing  routine  for  sparse  data  in  which 
directional  infonnation  was  retained  in  the  autocorrelation  function,  which  preserved  the 
2D  wavenumber  in  the  spectrum. 

B.  ICE  MOTION  DATA  PROCESSING 

Numerous  IABP  buoys  passed  through  the  area  of  the  BGEP  moorings  from 
2003-2005  (Fig.  11).  The  associated  raw  buoy  position  data  were  obtained  from  the 
IABP  office.  Each  file  represented  one  buoy  and  covered  a  period  of  about  1  year.  The 
following  processing  routine  was  applied  to  each  data  file  using  an  algorithm  written  in 
MATLAB. 
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Figure  1 1 .  Drift  tracks  of  all  buoys  in  mooring  area  from  2003-2005. 

First,  the  IABP  buoy  data  file  was  loaded  and  all  times  were  converted  to  year  day 
referenced  to  01  January  2003  as  the  start  point.  Next,  redundant  records  were  removed. 
The  record  was  then  scanned  to  find  times  when  neither  the  latitude  nor  longitude 
changed  for  successive  times.  These  records  were  removed  as  well.  Position  data  were 
then  converted  from  geographic  coordinates  of  latitude  and  longitude  to  a  polar 
stereographic  grid  compatible  with  the  International  Bathymetric  Chart  of  the  Arctic 
Ocean  (IBCAO).  This  coordinate  system  is  centered  with  the  origin  at  the  geographic 
North  Pole.  The  x-axis  measures  distance  in  meters  from  the  North  Pole  along  the  line  of 
90  degrees  East  longitude  and  the  y-axis  measures  distance  in  meters  from  the  North  Pole 
along  the  line  of  180  degrees  longitude.  The  buoys  did  not  report  position  data  on  a 
standard  time  schedule.  Therefore,  the  data  were  interpolated  onto  a  common  hourly 
time  vector  from  the  start  of  2003  to  the  end  of  2005.  This  allowed  buoy  positions  to  be 
available  for  every  hour  (on  the  hour).  The  polar  stereographic  buoy  positions  were 
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despiked  to  remove  anomalous  position  data.  A  remaining  problem  was  occasional,  long 
temporal  gaps.  Interpolating  across  these  would  have  introduced  erroneous  velocity 
estimates.  These  gaps  in  the  position  time  series  could  be  due  to  rime  building  up  on  a 
buoy’s  GPS  or  broadcast  antenna  at  some  point  and  then  getting  washed  off  later. 
Significant  data  gaps  were  identified  by  finding  position  fixes  that  occurred  more  than  6 
hours  apart.  Interpolated  data  between  these  fixes  were  discarded.  Velocity  components 

(u,v)  were  then  calculated  as  a  finite  difference  using  u=  —  =  ^ — — .  The  speed 

At  t2-tx 

(magnitude  of  the  overall  velocity)  was  determined  via  the  Pythagorean  theorem. 


The  next  task  was  to  use  the  processed  buoy  data  to  generate  a  time  series  of  u 
(east-west  velocities)  and  v  (north-south  velocities)  at  each  mooring.  The  following 
process  was  carried  out  for  each  mooring  location. 


First,  the  mooring’s  geographic  coordinates  of  latitude  and  longitude  were 
converted  to  the  polar  stereographic  grid  previously  described.  For  each  time  step  (every 
hour),  the  distance  between  all  contemporaneous  buoy  positions  and  the  mooring  were 
calculated.  All  buoys  within  250km  of  the  mooring  (Fig.  12)  at  that  time  were  flagged 
and  the  mean  of  their  velocity  components  (u,  v)  was  then  recorded  as  a  value 
representative  for  that  mooring.  A  time  series  was  created  at  each  mooring  by  repeating 
the  process:  advancing  one  time  step  at  a  time,  detennining  the  buoys  within  250km,  and 
recording  the  mean  of  their  motion  data. 
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Figure  12.  Tracks  of  all  IABP  buoys  passing  within  250km  of  mooring  C  from  2003- 
2005. 
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Ice  Motion  Time  Series 

Mooring  C:  August  2004-August2005 
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Figure  13.  Example  time  series  of  motion  compiled  at  Mooring  C:  August  2004- 
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Note  that  the  actual  direction  of  ice  motion  for  mooring  C  was  constantly 
changing  because  of  the  inertial  oscillations  (Fig.  13).  Therefore,  the  motion  data  needed 
to  be  segmented  into  a  number  of  segments  so  that  each  piece  was  approximately 
unidirectional.  For  the  entire  time  series  of  motion  (at  each  mooring),  the  direction  was 

f  v 

calculated  from  u  and  v  using  a  four  quadrant  arctangent  function  6  =  tan  —  .  Hourly 

\u) 

change  in  direction  as  well  as  cumulative  change  in  direction  were  derived.  Also,  the 
hourly  distance  travelled  and  the  cumulative  distance  travelled  were  calculated.  If  the 
cumulative  change  in  direction  or  the  cumulative  distance  travelled  exceeded  a  threshold 
(defined  here  as  45  degrees  for  direction  and  2,000m  for  distance),  the  segment  would 
end  with  the  point  previous  in  time,  so  as  not  to  exceed  the  threshold.  The  cumulative 
counters  were  reset  to  zero  and  the  process  was  repeated.  Segment  start  and  stop  times, 


and  average  speed  and  direction  were  recorded. 
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A  limitation  of  this  approach  is  that  it  assumes  that  the  ice  motion  is  continuous  in 
space  between  the  buoys  and  moorings,  ignoring  spatial  gradients  of  velocity.  For 
combining  mooring  and  motion  data,  Melling  et  al  (1995)  cautioned  against  using 
Lagrangian  methods  such  as  radio  beacons  or  GPS  receivers  in  buoys  mounted  in  the  ice 
because  they  quickly  leave  the  moorings  behind  and  do  not  give  true  motion  at  the 
mooring.  The  error  is  particularly  large  near  coastal  boundaries  and  during  the  ice  melt 
period.  In  this  case,  well  away  from  land  boundaries,  it  is  assumed  that  the  primary 
motion  is  due  to  inertial  oscillations  resulting  from  synoptic  scale  meteorological  events 
of  1,000-2,500  km.  With  forcing  of  this  scale,  it  is  assumed  that  ice  motion  is  continuous 
for  a  distance  of  250  km. 

C.  COMBINING  THE  ICE  DRAFT  AND  MOTION  DATA 

Without  motion  data,  all  that  is  recorded  at  the  ULS  is  a  time  series  of  ice  draft. 
The  addition  of  motion  data  helps  to  transfonn  this  time  series  into  a  2D  spatially 
referenced  description  of  ice  draft.  The  process  used  below  was  first  outlined  by  Hibler 
(1972). 

For  each  mooring,  the  ice  draft  information  was  loaded  and  corresponding  times 
were  converted  to  year  days  referenced  to  01  January  2003  to  be  consistent  with  previous 
data.  Using  each  time  segment  derived  from  the  motion  data,  the  corresponding  ice  draft 
time  series  had  its  mean  removed  and  the  autocorrelation  function  was  calculated  using 

1  N- 1 

Rxx  (Ax)  =  —£/(*/ )  fix,  +  At)  ,  where  Rxx  is  the  autocorrelation  value  (in  units  of  m 2 ), 

N  i= o 

N  is  the  total  number  of  data  points,  fix.)  is  the  value  of  the  draft  function  (in  units  of 
in  units  of  m)  evaluated  at  x  (in  units  of  m),  and  Ax  is  the  spatial  lag  (in  units  of  m)  (Fig. 
14). 
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Figure  14.  Ice  draft  time  series,  distance  series,  and  autocorrelation  function  for 
Segment  #50  at  Mooring  C. 


The  spatial  lags  Ax  were  a  function  of  segment  velocity.  These  lags  were 
calculated  using  Ax  =  At-  \lu2  +v2  ,  where  At  was  the  sampling  interval  (2s),  and 
\lu2  +  v2  was  the  segment  average  speed  (m/s). 

For  each  mooring,  all  segment  derived  autocorrelation  functions  were  then  plotted 
together  as  functions  of  lag  and  segment  direction.  The  result  was  a  series  of 
autocorrelation  functions  emanating  radially  from  the  mooring  like  spokes  of  a  wheel 
(Fig.  15). 
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2-D  Autocorrelation  Functions 

Mooring  C:  Aug  2004-Aug  2005 


Lag  in  X  Direction  (m) 


Figure  15.  Combined  autocorrelation  spatial  coverage. 

It  is  important  to  note  that  the  process  used  here  departed  from  the  one  outlined 
by  Hibler.  Hibler  recommended  normalizing  the  autocorrelation  functions.  However, 
normalization  strips  away  the  total  variance,  which  is  critical  for  interpreting  spectra. 
Normalized  autocorrelation  functions  are  unable  to  tell  the  difference  between  small  and 
large  roughness  elements.  The  large  elements,  the  ice  keels,  need  to  be  distinguished. 
Their  effects  are  distinct  from  roughness  elements  of  a  much  smaller  centimeter  scale.  A 
field  of  small  roughness  elements  oriented  in  one  direction  could  dominate  a  normalized 
2D  spectrum,  even  though  there  was  a  more  sparse  field  of  very  large  roughness  elements 
at  a  different  orientation  which  in  reality  dominates  the  physical  processes.  Hence,  in  this 
study,  variance  was  conserved  (i.e.  the  autocorrelations  were  not  normalized)  in  order  to 
make  direct  comparisons  between  different  locations,  times,  and  directions. 

It  is  also  important  to  note  that  this  step  of  combining  the  autocorrelation 
functions  is  crucial  because  it  mitigates  the  effects  of  sampling  across  features  at  an 
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angle.  The  relative  change  in  autocorrelation  functions  between  high  and  low  lag 
numbers  is  detennined  by  the  direction  that  the  linear  ridge  feature  is  crossed.  By  using 
the  autocorrelation  function  and  referencing  it  with  its  direction,  the  subsequent  2D 
spectrum  does  not  become  reddened  because  the  2D  wavenumber  is  conserved.  Energy 
peaks  will  occur  at  the  correct  wavenumbers  even  though  the  field  is  sparsely  sampled 
This  is  demonstrated  with  test  cases  in  section  III  A. 

This  combined  autocorrelation  function  is  an  intennediate  step.  The  desired  final 
product  is  the  2D  spectrum.  Since  the  autocorrelation  function  and  the  spectrum  are  a 
Fourier  transform  pair,  it  is  straightforward  to  go  from  one  to  the  other  using  the  Fast 
Fourier  Transform  (FFT).  However,  the  data  needs  to  be  on  a  continuous  grid  in  order  to 
use  the  FFT.  The  autocorrelation  function  had  some  gaps  in  coverage  and  so  needed  to 
be  processed  to  be  ready  for  the  2D  FFT. 

The  next  step  was  to  combine  all  of  the  autocorrelation  functions  by  binning 
them  in  polar  coordinates.  The  bins  chosen  were  10  degrees  wide  and  5m  long  to  achieve 
satisfactory  point  density.  It  became  apparent  that  each  autocorrelation  function  had  to 
be  interpolated  onto  a  common  vector  in  order  to  keep  the  spacing  of  the  points  constant 
when  compared  to  the  other  functions.  It  was  observed  that  due  to  the  fixed  sonar 
sampling  rate  of  2s,  autocorrelation  functions  corresponding  to  periods  of  slow  ice 
motion  had  their  data  points  spaced  much  more  closely  together  than  did  functions  from 
periods  of  faster  ice  motion.  The  result  was  that  a  particular  bin  might  have  10  points 
from  the  slower  autocorrelation  function  and  only  1  point  from  the  faster  one,  even 
though  they  both  crossed  the  bin  at  the  same  angle.  The  value  subsequently  assigned  to 
the  bin  was  biased  in  favor  of  the  slower  moving  segment  of  data.  The  autocorrelation 
data  was  therefore  interpolated  onto  vectors  of  common  length  and  point  spacing, 
specifically  2000m  and  0.5m  respectively.  This  issue  of  oversampling  during  periods  of 
slow  ice  motion  and  undersampling  during  periods  of  fast  ice  motion  was  addressed  by 
Melling  and  Riedel  (1995)  who  showed  that  resampling  at  a  constant  spatial  rate  was 
more  appropriate  than  using  the  time  series  when  trying  to  derive  independent  data  for 
statistical  comparisons.  The  mean  was  taken  for  each  bin  and  its  value  was  assigned  to  a 
point  at  the  bin’s  center.  These  binned  data  were  then  interpolated  and  transferred  from 
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the  polar  coordinate  system  to  a  Cartesian  grid  measuring  from  -1300  to  1300  m  in  both  x 
and  y  directions  with  grid  points  every  5m.  The  2D  FFT  of  the  2D  autocorrelation 
function  was  taken  to  obtain  a  2D  spectra  using 

i  N- 1 JV— 1 

S(k,l)  =  ,A Vj)e  2m{k-‘Ax<+l'Ayj'>  ^  where  k  and  /  are  the  horizontal  and  vertical 

N  i= o  j= o 

wavenumbers,  Axr  and  Av .  are  the  lags  in  the  x  and  y  directions,  and  R  is  the  value  of  the 
binned  and  interpolated  2D  autocorrelation  function  described  above. 


The  highest  resolvable  wavenumber,  the  Nyquist  wavenumber  was 

kN  ist  = — - — ,  where  Ax  =grid  spacing  (5m).  This  corresponded  to  a  wavelength  of 
2- Ax 

1 0m,  which  was  expected  to  be  suffucient  to  characterize  the  rapidly  changing  portion  of 
the  under  ice  environment. 


Synthetic  data  was  created  from  a  known  field  in  order  to  test  all  phases  of  the 
process  (see  section  III  for  detailed  results).  The  process  showed  itself  much  more 
sensitive  to  variations  in  angular  rather  than  radial  bin  width.  Decreasing  angular 
resolution  caused  the  spectra  to  get  smeared  out.  Increasing  it  caused  critical  decreases  in 
point  density  in  the  inital  binning  process  due  to  the  sparse  nature  of  the  data.  Radial 
variations  had  little  effect  on  the  spectra.  Fiowever,  as  the  radial  resolution  increased, 
there  was  not  enough  point  density  to  fill  every  bin.  If  it  decreased,  the  ultimate 
resolution  of  features  occuring  with  higher  wave  numbers  was  decreased  and  roughness 
on  these  scales  cannot  be  determined.  In  the  end,  both  angular  and  radial  resolutions 
used  were  a  compromise  in  the  fact  that  they  had  to  remain  coarse  enough  to  deal  with 
the  sparse  data. 

The  size  of  the  2D  autocorrelation  area  taken  for  the  2D  FFT  was  detennined  by 
the  length  scale  of  the  phenomenon  occuring  with  the  lowest  wave  number.  This 
phenomenon  was  assumed  to  be  distinct  ice  keels.  Historical  data  shows  the  spacing 
between  distinct  ice  keels  has  over  a  99%  probability  of  being  less  than  1000m  (Fig.  16) 
(Wadhams  and  Davy,  1986).  Therefore,  the  output  from  the  FFT  would  need  to  be  able  to 

resolve  a  wave  number  of  0.001m'1.  Using  resolution  =  — - — ,  it  follows  that  for  a  field 

N- Ax 
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with  grid  spacing  Ax  equal  to  5m,  the  required  resolution  of  0.001m"1  can  only  be 
achieved  if  the  number  of  points  N  is  greater  than  or  equal  to  200.  This  forces  the  size  fo 
the  field  to  go  into  the  2D  FFT  to  have  a  radius  of  at  least  1000m  and  the  field  to  measure 
at  least  2000m  across. 


Figure  16.  Cumulative  distribution  functions  of  ice  keel  spacing.  Plot  (a)  is  derived 
from  narrow  beam  sonar  in  the  Beaufort  Sea:  open  circles  have  a  threshold  of  0m  and 
closed  circles  have  a  threshold  of  3m.  Plot  (b)  is  derived  from  wide  beam  sonar  in  the 
Eurasian  Basin:  open  circles  have  a  threshold  of  0m  and  closed  circles  have  a  threshold 
of  15m  (after  Wadhams  and  Davy,  1986). 
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Even  when  using  the  synthetic  data  of  idealized  cases,  there  were  artifacts  in  the 
spectra  from  the  processing  routine.  To  attempt  to  decrease  these  artifacts,  zero  mean 
random  noise  was  added  to  the  2D  autocorrelation  surface  before  the  2D  FFT  converted 
the  field  into  spectra.  The  noise  level  used  was  1%  of  the  standard  deviation  of  the  entire 
2D  autocorrelation  surface.  Even  with  this  improvement,  there  continued  to  be  noise  at 
very  low  wave  numbers  corresponding  to  long  wavelengths.  2D  convolution  filtering 
failed  to  make  the  plots  easier  to  interpret.  Much  of  the  noise  was  at  wave  numbers 
below  0.001m'1,  which  had  previously  been  set  as  the  maximum  length  scale  imposed  by 
independent  keel  spacing.  Therefore,  although  it  was  apparent  that  the  process  could  not 
differentiate  between  the  processing  artifacts  and  real  data  at  very  low  wave  numbers 
(long  wavelengths),  it  was  considered  to  be  acceptable  because  energy  is  not  expected  to 
be  there  due  to  the  physical  scaling  argument. 
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IV.  ANALYSIS  AND  DISCUSSION 


A.  TEST  CASES 

To  confirm  that  the  technique  worked  as  expected  and  to  aid  the  interpretation  of 
real  world  data,  a  number  of  test  cases  were  conducted  where  a  known  field  of  under  ice 
features  was  put  through  the  processing  routine.  The  random  fields  were  generated  and 
then  sampled  using  real  ice  motion  observations.  The  segments  used  were  from  Mooring 
C  from  August  2004-August  2005.  For  consistency,  the  Mooring  C  trajectories  were 
used  to  generate  ice  draft  data  from  the  known  fields  for  all  of  the  test  cases. 

1.  Single  Sine  Wave 

The  first  test  case  was  that  of  a  field  created  from  a  single  sine  wave  with 
amplitude  5m,  wavelength  200m,  and  oriented  so  the  crests  form  a  line  running  vertically 
across  the  field  (Fig.  17(a)).  When  put  through  the  processing  routine,  the 
autocorrelation  function  showed  the  orientation  correctly  (Fig.  18),  and  the  resulting 

spectra  put  the  peak  energy  at  the  correct  wave  number  and  direction,  [k,/]  =  [^^,0] 

(Fig.  17(b)  and  (c)). 
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Field  from  Single  Sinusoid 
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Figure  17.  (a)  Field  generated  from  single  sinusoid,  (b)  Resulting  2D  spectrum  with 

energy  peaks  at  correct  wavenumbers,  (c)  Spectral  cross  section  with  vertical  axis  scaled 
by  wavenumber.  Note  peak  energy  at  correct  wave  number. 
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2  D  Autoccrrelation  Surface 
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Figure  18.  2D  autocorrelation  plot  for  field  created  from  a  single  sinusoid.  Note  the 
distortion  caused  by  the  angular  binning  process  as  angles  come  less  perpendicular  to  the 
axis  of  the  crests. 


2.  Single  Sine  Wave  Derived  from  Different  Motion  Data 

The  second  test  case  was  another  single  sine  wave  with  identical  parameters  as  in 
Case  1.  The  object  of  this  case  was  to  determine  if  the  output  was  independent  of  the 
velocity  estimates  used  to  translate  the  ULS  time  series  into  a  spatial  series.  The  identical 
field  from  Case  1  was  generated  but  this  time  it  was  sampled  by  ice  motion  data  for 
Mooring  C  from  the  previous  year  (August  2003-August  2004).  Again,  the  spectra  put 
the  peak  energy  at  the  correct  wave  number,  suggesting  that  the  technique  is  insensitive 
to  the  details  of  the  conversion  of  the  time  series  to  a  spatial  series  with  the  buoy-derived 
velocities  (Fig.  19). 
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Field  from  Single  Sinusoid 
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Figure  19.  (a)  Field  generated  from  sinusoid,  (b)  2D  spectrum  from  field  sampled 

with  different  ice  motion  data,  (c)  Spectral  cross  section. 
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3. 


2  Sine  Waves  of  Different  Wavelengths  but  Different  Amps  and 
Different  Orientation 


This  test  case  sampled  a  field  made  by  combining  two  differing  sinusoidal  fields 
(Fig.  20).  The  first  field  had  a  wavelength  of  200m,  vertically  oriented  crests,  and  an 
amplitude  of  5m.  This  corresponded  to  a  wavenumber  of  k=.005  m"1  and  logl0(k)=-2.3. 
The  second  field  had  sinusoids  of  500m,  diagonally  oriented  crests,  and  amplitude  of 
20m.  This  corresponded  to  a  wavenumber  of  k=.002  m'1  and  logl0(k)=-2.69. 

At  a  glance,  it  may  appear  that  only  the  larger  amplitude  wave  was  represented  in 
the  2D  spectrum  (Fig.  21(b)).  However,  cross  sections  through  the  2  directions  of 
interest  show  that  both  sinusoids  were  correctly  resolved  (Fig.  21(c)).  Not  only  were 
peaks  located  in  the  correct  wave  numbers,  but  their  relative  sizes  were  consistent  with 
their  relative  amplitudes.  This  reflects  the  larger  variance  in  larger  amplitude  ice  keels. 
This  was  one  of  the  most  important  test  cases  because  it  showed  that  the  process  is  robust 
enough  to  properly  display  a  field  created  from  sinusoids  with  differing  orientation, 
wavelength,  and  amplitude. 


Sinusoid  #1 


Figure  20.  Sinusoids  #1  and  #2  which  combine  to  make  up  the  composite  field  of  two 
sinusoids  with  differing  amplitude,  wavelength,  and  orientation. 
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Composite  Field  from  Two  Sinusoids 


(b) 


2-D  Spectra  (square  root) 
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Figure  21.  (a)  Composite  field  from  two  sinusoids,  (b)  2D  spectrum  of  field,  (c) 

Spectral  cross  sections  with  peak  energy  and  relative  amplitude  correctly  represented. 
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4.  Single  Sine  Function 


This  test  case  was  to  test  how  well  the  process  worked  on  a  non-sinusoidal  field. 
A  field  was  created  with  a  sine  function  repeating  periodically  with  a  wavelength  of 
100m  (Fig.  22(a)).  To  approximate  the  shape  of  a  true  ice  keel,  the  amplitude  was  set  at 
5m,  it  was  20m  wide  base.  The  peak  amplitude  of  the  spectra  was  at  the  correct 
wavenumber  (Fig.  22(b)).  However,  hannonics  were  observed  in  the  spectra  (Fig.  22(c)). 
This  is  due  to  the  nature  of  the  2D  FFT  process,  which  struggles  to  fit  high  wavenumber 
sinusoids  into  the  observed  profile.  The  harmonics  were  expected  however,  as  the 
Fourier  transform  of  the  sine  function  is  the  square  wave.  The  implication  is  that  the 
process  can  correctly  represent  non-sinusoids. 
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Single  Sine  Function 
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(b) 


2-D  Spectra  (square  root) 
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Figure  22.  (a)  Field  from  single  sine  function,  (b)  2D  spectrum  from  field,  (c) 

Spectral  cross  section;  note  the  decreasing  harmonics  at  higher  wavenumbers. 
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B.  BEAUFORT  GYRE  DATA 

1.  Mooring  C:  August  2004-August  2005 

Mooring  C  had  the  most  dense  buoy  coverage  and  is  used  to  illustrate  the 
Beaufort  Gyre  data.  The  case  described  below  is  from  August  2004-  August  2005,  the 
second  year  of  deployment.  This  dataset  had  the  best  coverage  with  respect  to  the  2D 
autocorrelation  function  and  therefore  is  considered  the  most  robust. 

The  2D  spectrum  (presented  as  the  square  root  of  variance  for  display  purposes) 
appears  by  eye  to  be  directionally  variable  (Fig.  23a).  Eight  spectral  cross  sections  taken 
in  directions  ranging  from  000  to  157.5°  in  22.5°  intervals  show  that  the  spectra  have  a 
power  law  dependence  on  radius  (Fig.  23(b)  and  (c)).  Most  of  the  variance  of  the 
underside  morphology  is  contained  in  a  wavenumber  band  between  0.001  and  0.01m"1 
(100- 1000m).  Spectral  levels  ranged  from  1.2  to  4.2  x  104  m2.  There  is  no  distinct, 
preferred  wavelength. 
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Figure  23.  Mooring  C:  2004-2005.  (a)  2D  spectrum,  (b)  8  energy  conserving 

spectral  cross  sections  (000  to  157.5°  each  22.5°).  (c)  8  spectral  cross  sections. 
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In  order  to  determine  whether  or  not  the  cross  sections  are  statistically  different 
from  one  another,  the  standard  deviation  needed  to  be  determined.  This  allowed 
confidence  limits  to  be  plotted  against  one  another  to  check  for  overlapping  regions, 
which  would  mean  that  there  is  not  a  difference.  However,  if  the  regions  bounded  by  the 
confidence  intervals  were  mutually  exclusive,  they  are  different  and  indicate  anisotropy. 

Determination  of  the  confidence  intervals  was  not  straightforward.  Since  the  data 
was  taken  at  one  frequency,  autocorrelated,  subsampled,  binned,  meaned  and  gridded  in 
autocorrelation  space  before  passing  through  a  2D  FFT,  degrees  of  freedom  are 
impossible  to  determine  analytically.  Therefore,  the  degrees  of  freedom  were  detennined 
empirically  by  Monte  Carlo  simulation.  A  field  was  created  which  would  ultimately 
approximate  the  characteristics  of  the  observations.  This  was  done  by  starting  with  a 
field  of  random  numbers,  taking  the  2D  FFT,  multiplying  by  1/k  and  then  taking  the 
inverse  2D  FFT,  to  get  back  to  a  Cartesian  field.  This  new  field  had  the  properties  that 
when  it  was  put  through  the  process  it  yielded  a  2D  spectrum  which  was  comparably  red 
like  the  real  data,  but  composed  of  random  numbers. 

This  random  field  was  sampled  using  real  buoy  motion  and  put  through  the 
processing  routine.  A  radial  cross  section  was  taken  through  the  resulting  spectrum.  This 
entire  process,  from  creating  the  random  field  through  taking  the  cross  section  was 
repeated  30  times.  For  each  wavenumber  bin,  the  standard  deviation  of  the  cross  sections 
was  determined.  A  line  was  fit  to  approximate  the  standard  deviation  as  a  function  of 
wavenumber.  This  represented  the  sigma  value  used  for  confidence  limits. 

Two  cross  sections  were  taken  in  directions  which  appeared  to  be  similar  (Fig. 
24(a)).  Even  at  1  sigma,  the  spectra  showed  no  regions  of  mutual  exclusivity  (Fig. 
24(b)).  This  suggested  isotropy.  At  2  sigma,  the  confidence  intervals  are  even  wider  and 
the  overlap  even  greater  (Fig.  24(c)).  The  result  was  to  conclude  that  the  spectrum  is 
isotropic  in  these  two  directions. 

Next,  2  cross  sections  were  taken  in  directions  which  appeared  to  be  significantly 
different  (Fig.  25(a)).  At  1  sigma,  the  spectra  did  show  regions  of  mutual  exclusivity, 
especially  at  the  high  wavenumbers,  suggesting  anisotropy  (Fig.  25(b)).  However,  these 
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regions  of  exclusivity  disappeared  at  all  except  for  a  few  points  when  the  confidence 
intervals  were  set  to  2  sigma,  which  suggests  isotropy  (Fig.  25(c)).  The  result  was  to 
conclude  that  there  is  ‘weak’  anisotropy  in  this  spectrum.  This  is  defined  to  mean  that 
the  results  were  isotropic  except  for  low  levels  of  statistical  confidence. 
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Figure  24.  (a)  2  radial  slices  through  2D  spectrum,  (b)  1  sigma  confidence  intervals 

for  each  slice.  The  red  lines  represent  the  +  and  -  1  sigma  confidence  intervals  for  the 
spectral  slice  taken  along  the  red  dotted  line  in  (a)  and  the  green  lines  are  the  confidence 
intervals  for  the  slice  taken  along  the  green  dotted  line,  (c)  2  sigma  confidence  intervals. 
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Figure  25.  (a)  2  radial  slices  through  2D  spectrum,  (b)  1  sigma  confidence  limits  for 

each  slice.  The  red  lines  represent  the  +  and  -  1  sigma  confidence  intervals  for  the 
spectral  slice  taken  along  the  red  dotted  line  in  (a)  and  the  green  lines  are  the  confidence 
intervals  for  the  slice  taken  along  the  green  dotted  line,  (c)  2  sigma  confidence  intervals. 
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To  ensure  that  the  observed  anisotropy  was  not  a  function  of  the  binning  process, 
a  “red”  field  like  the  one  used  to  generate  the  confidence  intervals  was  sent  through  the 
process.  Cross  sections  were  taken  in  the  2  directions  which  showed  the  greatest 
amounts  of  anisotropy  for  the  previous  run.  The  result  showed  that  the  spectrum  was 
isotropic,  even  at  one  sigma  (Fig.  26).  This  shows  that  the  observed  isotropy  was  a 
characteristic  of  the  real  data  and  not  the  processing  routine. 
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Figure  26.  1  sigma  confidence  intervals  from  the  spectral  cross  sections  from  the 

synthetic  “red”  field  taken  in  the  same  directions  as  Fig.  25. 

2.  Summary  for  Other  Data  Sets 

There  were  directional  differences  observed  in  the  spectra  of  all  data  sets  at  1 
sigma.  These  suggest  anisotropy.  However,  this  was  never  true  at  2  sigma,  where  the 
spectra  appeared  isotropic.  Furthermore,  although  the  isotropic  regions  observed  usually 
occurred  at  high  wavenumbers,  they  were  also  sometimes  observed  throughout  the  rest  of 
the  spectra.  As  for  the  data  from  Mooring  C  from  2004-2005,  most  of  the  variance  of  the 
other  datasets  is  contained  in  a  wavenumber  band  between  0.001  and  0.01m'1  ( 1 00- 
1000m).  Again,  there  are  no  distinct,  preferred  wavelengths.  Spectral  levels  ranged  from 
0.2  -  1.6  x  104  m2  for  Mooring  A  from  2004-2005  to  0.3  -  3.3  x  104  nr  for  Mooring  C 
from  2003-2004. 
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In  summary,  the  spectra  are  weakly  anisotropic  for  the  wavenumbers  resolved  in 
this  study. 

C.  SPATIAL  AND  TEMPORAL  COMPARISONS 

The  weak  anisotropy  observed  was  not  consistent  from  one  mooring  to  another. 
The  2D  spectra  are  displayed  in  Figs.  27  and  28.  Spatial  and  temporal  comparisons  failed 
to  detennine  a  consistent  directional  preference.  However,  most  of  the  variance  was 
contained  in  the  .001  to  .01m'1  wavenumber  band  (lOO-lOOOm)  for  all  datasets.  Spectral 
levels  varied  significantly  between  the  datasets. 
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Figure  27. 
Mooring  A. 
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2D  spectra  of  under  ice  roughness  from  August  2003-August  2004. 
(b)  Mooring  B.  (c)  Mooring  C. 
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D. 


COMPARISONS  TO  OTHER  RESEARCHERS 


This  weak  anisotropy  does  not  agree  with  the  findings  of  Davis  and  Wadhams 
(1995).  While  working  from  a  submarine,  they  gathered  concurrent  upward  looking 
echosounder  and  towed  sidescan  sonar  which  allowed  for  determination  of  both  ice  draft 
and  ice  keel  orientation.  Their  findings  showed  that  there  was  anisotropy  in  the 
alignment  of  the  large  ice  roughness  elements,  the  ice  keels.  Their  data  showed  lineation 
which  fell  roughly  into  two  preferred  alignments.  When  considered  in  the  context  of  the 
area  in  which  the  data  were  gathered,  the  directions  are  roughly  along  and  perpendicular 
to  the  Greenland  coast  (Fig.  29).  The  observed  ice  keels  could  be  the  result  of  ridging 
due  to  compression  caused  by  downstream  boundary  and  shearing  caused  by  motion  of 
Transpolar  Drift  Stream. 


Figure  29.  Ridge  orientation  vectors  (Davis  and  Wadhams,  1995). 

Comparisons  to  the  findings  of  Davis  and  Wadhams  and  this  project  should  be 
made  carefully.  This  project  used  data  from  the  Beaufort  Gyre,  well  away  from  land 
boundaries.  Furthermore,  the  data  was  collected  from  moorings  in  deep  water  (greater 
than  3500m)  so  there  was  no  bottom  interaction  or  grounding  to  cause  a  preferred  angle 
of  orientation  for  roughness  elements. 
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Hibler  (1972)  first  suggested  the  process  upon  which  this  project  was  based. 
Using  his  version  of  data  processing,  he  found  anisotropy  to  suggest  a  preferred  lineation 
of  pressure  ridge  sails  in  the  Beaufort.  However,  this  was  based  on  one  case  using  three 
intersecting  laser  profiles  about  2400m  long  collected  by  aircraft.  In  contrast,  this  project 
used  several  hundred  profiles  in  three  locations  in  the  Beaufort  over  the  course  of  two 
years.  Again,  comparisons  must  be  made  carefully. 

E.  TURBULENCE  IMPLICATIONS 

The  implications  of  the  weak  anisotropy  that  was  observed  is  that  directionality  is 
important  at  very  small  scales.  In  turbulence  simulations,  where  roughness  is  significant, 
a  ‘spectral’  approach  is  most  appropriate  to  approximate  a  near-isotropic  field.  In 
contrast,  most  conceptual  models  of  flow  over  roughness  have  been  developed  with 
regular  arrays  of  ridges  and  keels  with  a  single  orientation  and  fixed  spacing  (Arya, 
1975). 

Density  profiles  obtained  at  APLIS-2007  in  the  Beaufort  Sea  showed  pycnoclines 
which  were  significantly  shallower  than  those  historically  observed  (Fig.  30).  This 
shoaling  of  the  thennocline  is  possibly  explained  by  the  recent  accelerating  rate  of  ice 
melt.  Pycnoclines  at  shallow  depths  can  produce  the  “dead  water”  effect  if  an  ice  keel  is 
dragged  through  them.  With  any  anisotropy  in  the  ice  roughness,  especially  for  deep  ice 
keels  whose  depth  is  comparable  to  or  exceeding  the  pycnocline  depth,  there  will  be 
significant  effects  on  the  momentum  budget.  Energy  will  be  lost  at  variable  rates 
depending  on  the  direction  that  the  ice  flow  is  drifting. 


50 


Yearday  82.37 

Yearday  8-4.42 

Yearday  87.66 

16 

Id- 

i: 

20 

20 

•mV  \ 

30 

Id- 

mV  \  -I 

40 

40 

•40  r 

1  „ 

M 

M 

J 

v> 

\  1 

-*[  \  4 

70 

70*-  \  4 

70‘  V  4 

00 

Nr  \  j 

o:  i 

>0 

00-  H 

00r 

100 

1422  1023  1024  1*2&  1026 

M  1022  1023  1024  1029  102« 

00  1022  1023  1024  10»  102 

Density  (kg/m*) 

Density  (kg/m1) 

Density  (kg/nf) 

Figure  30.  Density  profiles  from  APLIS  2007.  Note  the  shallow  pycnocline; 
historically  the  pycnocline  has  been  observed  near  35m. 

F.  POSSIBLE  MITIGATING  EFFECT 

An  observation  made  recently  showed  that  the  ice  floe  used  for  APLIS-2007 
rotated  significantly,  changing  heading  by  11  degrees  in  7  days  (Fig.  31)  (Hutchings, 
2007).  This  has  the  potential  to  mitigate  the  effects  of  weak  directional  dependence  that 
were  observed  in  this  study,  making  the  time  averaged  roughness  isotropic.  The  APLIS- 
2007  observation  is  considered  abnormal  when  compared  to  previous  ice  camp  rotation 
measurements.  However,  as  the  Arctic  pack  ice  becomes  thinner  and  more  mobile,  more 
rotation  is  possible.  This  would  make  conditions  more  isotropic,  and  should  be 
investigated  further  to  determine  if  it  is  a  new  trend. 
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Figure  31. 


APLIS  Rotation 
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Plot  of  APLIS  2007  ice  floe  rotation  (Hutchings,  2007). 
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V.  CONCLUSIONS 


The  roughness  of  the  ice  underside  plays  a  role  in  turbulent  ocean  heat,  salt  and 
momentum  fluxes  that  determine  the  delicate  balance  between  surface  forcing  and  the 
ocean  interior  that  sustains  or  reduces  the  perennial  ice  cover.  It  has  climate, 
engineering,  and  naval  implications  as  well. 

This  project  took  sparsely  sampled  ice  draft  data  and  processed  it  in  such  a  way 
that  the  2D  wavenumbers  were  conserved  in  order  to  produce  a  realistic  2D  spectrum  of 
under-ice  roughness.  The  process  was  validated  using  several  test  cases.  It  showed  weak 
anisotropy  of  under-ice  roughness  at  scales  between  10m  and  2000m  for  pack  ice  in  the 
Beaufort  Gyre.  Error  estimates  were  determined  by  Monte  Carlo  simulations.  For  each 
dataset,  most  of  the  variance  of  the  underside  morphology  was  contained  in  a 
wavenumber  band  between  0.001  and  0.01m"1  (lOO-lOOOm).  There  were  no  distinct, 
preferred  wavelengths.  Spectral  levels  ranged  from  1.2  -  4.2  x  104  m2  for  Mooring  C 
(2004-2005)  to  0.2  -  1.6  x  104  m2  for  Mooring  A  (2004-2005). 

Future  studies  will  benefit  from  the  new  mooring  configurations  in  the  Beaufort 
Gyre  Exploration  Project.  The  new  moorings  have  colocated  UFS  and  ADCPs  which 
when  set  in  bottom  tracking  mode  give  the  true  local  ice  motion  at  the  mooring.  There  no 
longer  needs  to  be  interpolation  between  buoys  or  assumptions  of  continuous  ice  motion. 
This  increased  amount  of  data  will  allow  tighter  limits  in  the  processing  routine  and 
ultimately  a  higher  resolution  2D  roughness  spectrum. 
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